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Abstract 

We present a novel way of constructing reduced models for sys- 
tems of ordinary differential equations. In particular, the approach 
combines the concepts of renormalization and effective field theory de- 
veloped in the context of high energy physics and the Mori-Zwanzig 
formalism of irreversible statistical mechanics. The reduced models 
we construct depend on coefficients which measure the importance of 
the different terms appearing in the model and need to be estimated. 
The proposed approach allows the estimation of these coefficients on 
the fly by enforcing the equality of integral quantities of the solution 
as computed from the original system and the reduced model. In this 
way we are able to construct stable reduced models of higher order 
than was previously possible. The method is applied to the problem of 
computing reduced models for ordinary differential equation systems 
resulting from Fourier expansions of singular (or near-singular) time- 
dependent partial differential equations. Results for the ID Burgers 
and the 3D incompressible Euler equations are used to illustrate the 
construction. We also present, for the ID Burgers and the 3D Euler 
equations, a simple and efficient recursive algorithm for calculating the 
higher order terms. 



Introduction 



Spatial discretizations or Fourier expansions of the solutions of partial differ- 
ential equations (PDEs) which depend on time lead to systems of ordinary 
differential equations (ODEs). The most difficult case arises when the solu- 
tion of a PDE becomes singular in finite time. At such instants the solution 
of the PDE develops activity down to the zero length scale. A brute force 
numerical simulation (no matter how large) of such a solution is bound to 
fail because the simulation has a finite resolution and thus will be unable 
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to resolve all the length scales down to the zero scale. When the solution 
develops activity at a scale smaller than the smallest scale available to the 
simulation, the numerically computed solution becomes underresolved. This 
leads to a rapid deterioration of the accuracy of the simulation. 

The notion of propagation of activity to smaller and smaller scales de- 
pends on the physical context of the PDE. In some cases, e.g. the 3D Euler 
or Navier-Stokes equations [1J, this could mean a cascade of energy to smaller 
and smaller scales. In other cases, e.g. the nonlinear focusing Schrodinger 
equation [2], this could mean a cascade of mass to smaller and smaller scales. 
Irrespective of the specific physical context, the problem facing the numeri- 
cal analyst is how to use a finite simulation and yet prevent the computed 
solution from suffering a loss of accuracy. In other words, how to construct 
a numerical method which reproduces correctly the features of the solution 
of the original equation at the length scales that are available numerically. 
This is the motivation behind the construction of reduced models (see e.g. 

mm)- 

By construction, a reduced model must allow for energy (mass) to escape 
from the scales that are accessible to the simulation (called resolved scales or 
modes) to the inaccessible scales (called unresolved). The main difficulty in 
constructing an accurate reduced model is the need to estimate the correct 
rate at which activity is propagated from the resolved to the unresolved 
scales. The Mori-Zwanzig (MZ) formalism [5j [6] proceeds by dividing the 
available resolution into resolved and unresolved parts. Then, it constructs 
a reduced model for the resolved scales and uses the unresolved scales to 
effect the drain of energy (mass) out of the resolved scales. 

Although the MZ formalism allows for the construction, in principle, of 
an exact reduced model it has two drawbacks (which are also shared by any 
other reduction formalism). First, the reduced model can be, in general, 
prohibitively expensive to calculate. The reason is that one must obtain an 
accurate representation of the behavior of the unresolved scales before they 
can be safely eliminated. Obtaining this representation can be rather costly. 

The second drawback is more subtle and has not been adequately ap- 
preciated by the scientific computing community. It is specific to the case of 
constructing reduced models for singular PDEs or in general for systems of 
ordinary differential equations which are larger than any available numerical 
resolution. Suppose that you have to construct a reduced model of a full 
system which is larger than any available numerical simulation. Let us call 
this system SI. Exactly because SI is larger than any available numerical 
simulation, if we want to construct a reduced model we have to use as a 
starting point a system, call it S2, whose size is smaller than the size of SI. 
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Suppose that you start with S2 and use the MZ formalism (or any other 
reduction formalism for that matter) and construct an exact reduced model 
S3 for a subset of S2. An exact reduced model means that if one evolves 
S2 and S3 separately, then the behavior of the scales resolved by S3 will be 
the same as the behavior of the scales in S3 predicted by the simulation of 
the system S2. However, and this is the heart of the problem, since S2 itself 
will become eventually underresolved, the exact reduced model S3 will also 
become underresolved. In other words, the predictions of the exact reduced 
model S3 can only be trusted for as long as the predictions of the system S2 
can be trusted. As a result, any reduced model that has any chance of being 
accurate for longer times cannot be exact. 

There are examples of inexact reduced models, e.g. the i-model [H El 
[U E], coming from the MZ formalism, which have been applied to singular 
PDEs and shown numerically to be relatively accurate for long time inter- 
vals. However, the t-model's accuracy is difficult to assess beforehand and 
the reason for its relative success has remained a mystery (see also [TO] for an 
application to the 3D Navier-Stokes equations which shows that the i-model, 
while not bad, is in need of some modification). In order to construct better 
reduced models we need to incorporate dynamic information from the full 
system which will help us decide which of the terms appearing in the exact 
reduced model are the ones that are most important. In this way, we can 
construct an inexact but accurate reduced model by keeping the important 
terms and disregarding the unimportant ones. The underlying principle of 
this approach is the same as the principle behind the construction of effective 
field theories in high energy physics [11] . 

The way we propose to address the problem of constructing better re- 
duced models is to embed the MZ reduced models in a larger class of reduced 
models which share the same functional form as the MZ reduced models but 
have different coefficients in front of the various terms that appear in the re- 
duced models. Then, one can estimate these coefficients on the fly while the 
original system of equations is still valid. The estimation of the coefficients 
is achieved by requiring that certain integral quantities (e.g. l p norms) in- 
volving only resolved scales, should acquire the same values when computed 
from the original system and the reduced model. The constraints used to 
obtain the coefficients are the analog of the "matching conditions" used in 
effective field theory. Also, the approach is the time-dependent analog of the 
process of renormalization used in high energy and condensed matter physics 
[12\ I13j . Before the original system ceases to be valid, one reverts to the 
reduced model with the various coefficients having their estimated values. 
We call the proposed approach the renormalized Mori-Zwanzig algorithm. 
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In addition to computing the coefficients of the reduced model on the fly, 
the proposed approach allows us to gauge the importance of the different 
terms appearing in the reduced model. While at the moment a rigorous proof 
is lacking, the numerical results (see Section 12. 3p suggest that the values of 
the coefficients depend on the ratio of the smallest scale present in the initial 
condition to the smallest scale resolved by the reduced model. The larger 
is this ratio, the larger are the coefficients for more terms and thus, more 
terms are needed in the reduced model. On the other hand, for the case 
when the ratio is small only a few number of terms need to be kept without 
sacrificing the accuracy of the model. As we have already said, this situation 
is analogous to the principle underlying the effective field theory method in 
high energy physics [11]. It is also interesting to investigate to what extent 
the values of the coefficients are determined by the scaling symmetries of 
the underlying PDE. This could allow us to identify different universality 
classes of PDEs based on their scaling symmetries (see also Section 12. 3ft . 

Finally, we note that the proposed method of computing the coefficients 
was first presented by the author in [T3]. The goal in [T5] was to construct 
a mesh refinement scheme to allow us to reach the singularity instant more 
efficiently. Then, the computed coefficients were used in a fixed point anal- 
ysis of the singularity in order to estimate the blow-up exponents. In the 
current work the goal is different. We not only want to reach the singularity 
instant or, more generally, the time that the full system runs out of resolu- 
tion, but also follow the solution for later times. Also, in the current work, 
we present a new way of expanding the memory term of the MZ formalism. 
This allows us, at least for Burgers and Euler, to calculate recursively and 
efficiently (and with minimal storage requirements) the higher order terms 
in the memory expansion. This is a welcome feature since analytical com- 
putation of the higher order terms quickly becomes exhausting and, thus, 
prone to errors. 

The paper is organized as follows. Section [1] presents the main construc- 
tion along with a brief presentation of the MZ formalism and how to compute 
the coefficients of the reduced model on the fly. Section [2] presents numer- 
ical results for the ID Burgers and the 3D incompressible Euler equations. 
Section [3] contains a discussion of the results and presents some directions 
for future work. 
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1 Renormalization of Mori-Zwanzig reduced mod- 
els 



In this section we explain the main ideas behind the proposed approach. In 
Section 11.11 we set up the notation for the original system and the reduced 
model in an abstract way which does not make reference to any specific 
method for obtaining the reduced model. In Section 11.21 we show how to 
obtain the coefficients for the reduced model. In Section 11.31 we give a brief 
presentation of the MZ formalism which allows us to obtain the functional 
form of the terms appearing in the reduced model. In Section fL^l we combine 
the ideas in Section 11.21 with the MZ formalism from Section 11.31 to derive 
the proposed algorithm for computing renormalized MZ reduced models. 



1.1 Full and reduced systems 

Suppose that we want to construct a reduced model for the partial differen- 
tial equation (PDE) 

v t + H(t,x,v,v x , ...) = 

where H is a, in general nonlinear, operator and x € D C R rf (the construc- 
tion extends readily to the case of systems of partial differential equations). 
After spatial discretization or expansion of the solution in series, the PDE 
transforms into a system of ordinary differential equations (ODEs). For 
simplicity we restrict ourselves to the case of periodic boundary conditions, 
so that a Fourier expansion of the solution leads to system of ODEs for the 
Fourier coefficients. To simulate the system for the Fourier coefficients we 
need to truncate at some point the Fourier expansion. Let FUG denote the 
set of Fourier modes retained in the series, where we have split the Fourier 
modes in two sets, F and G. We call the modes in F resolved and the modes 
in G unresolved. The reduced model involving only the resolved modes F 
will be called the reduced system and the system involving both the resolved 
and unresolved modes F U G will be called the full system. 

The main idea behind the algorithm is that the evolution of moments of 
the reduced set of modes, for example l p norms of the modes in F, should 
be the same whether computed from the full or the reduced system. This 
requirement will eventually allow us to compute the coefficients appearing 
in the reduced model (see Section 1 1 . 2f) . 

The full system of equations for the modes F U G is given by 

^ = *«,.<«)), 
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where u = ({uk}), k € F U G is the vector of Fourier coefficients of u 
and R is the Fourier transform of the operator H. The system should be 
supplemented with an initial condition u(0) = uq. The vector of Fourier 
coefficients can be written as u = (u,u), where u are the resolved modes 
(those in F) and u the unresolved ones (those in G). Similarly, for the 
right hand sides (RHS) we have R(t,u) = (R(t,u),R(t,u)). Note that the 
RHS of the resolved modes involves both resolved and unresolved modes. In 
anticipation of the construction of a reduced model we can rewrite the RHS 
as R{t,u) = R<P)(t,u) = {R^{t,u),R^{t,u)). 

In general, when one constructs a reduced model, additional terms ap- 
pear on the RHS of the equations of the reduced model (see Section 11.31 for 
more details). The role of these additional terms is to account for the in- 
teractions between the resolved and unresolved modes, since the unresolved 
modes no longer appear explicitly in the reduced model. As is standard in 
renormalization theory [15] . one can augment the RHS of the equations in 
the full system by including such additional terms. That is accomplished 
by multiplying each of these additional terms by a zero coefficient. In this 
way, the reduced and full systems' RHSs have the same functional form. In 
particular, for each mode Uk, k 6 F U G, we can rewrite R k °\t,u) as 

m 

Rf\t,u(t)) = Y, a ? )R ik( t > «(<)). 

i=l 

where R^(t,u(t)) = R k ° ] (t,u(t)) and R$(t,u(t)), for i = 2,...,m are of 
the same functional form as the additional terms which appear in the reduced 
model. This is easy to do by taking = 1 and =0, for % = 2, . . . , m. 
Thus, the equation for the the mode u&, k € F U G is written as 

du ^) _ o. n - ni°)u ,.r+\\ - „(°) p(°)/ 



R k (t,u) = R k °\t,u(t)) = Y,af )) R? k \t,u(t)) (1) 



dt 

i=i 

Correspondingly, the reduced model for the mode u' k , k € F is given by 

^=Ri 1 \t,At)) = ±a^R^(t,At)) (2) 

i=l 

with initial condition u' k (0) = iiQk- We repeat that the functions RjJ} , i = 
1, . . . , m, k € F, have the same form as the functions R^) , i = 1, . . . , m, k 6 
F, but they depend only on the reduced set of modes F. Dimensional reduc- 
tion transforms the vector = (af ,...,%') to = (a[ ,...,Om ). 
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This allows one to determine the relation of the full to the reduced system 
by focusing on the change of the vector to Also, the vectors and 
g^ 1 ** do not have to be constant in time. This does not change the analysis 
that follows. 

Define m quantities Ei, i = 1, . . . ,m involving only modes in F. For 
example, these could be l p norms of the reduced set of modes. To proceed 
we require that these quantities' rates of change are the same when computed 
from (P) and (J2J), i.e. 

dEi(u) dEi(u') . 

Note that similar conditions, albeit time-independent, lie at the heart of the 
renormalization group theory for equilibrium systems ([15] p. 154). Also, 
the conditions ([3]) are the analog of the "matching conditions" underlying 
the construction of effective field theories |llj . 

1.2 How to compute the coefficients of the reduced model 

When we only know the functional form of the terms appearing in the re- 
duced model but not their coefficients it is not possible to evolve a reduced 
system. We present a way of actually computing the coefficients of the re- 
duced model as needed. If the quantities Ei, i = 1, . . . , vn are e.g. lp norms 
of the Fourier modes, then we can multiply Equations (|2]) with appropriate 
quantities and combine with Equations (|3|) to get 

dEi(u) _ ^ (i) 



d,t 



dE 2 (u) _ (i) 



dt 



dE m (u) 



i=l 



m 



dt ^ 

i=i 



af> B mi {t,u{t)) 



da) 



where Bij = — I dEl ^ - ) , i,j = l,...,m are the new RHS functions that 



appear. Note that the RHS of the equations above does not involve primed 
quantities. The reason is that here the reduced quantities are computed by 
using the values of the resolved modes from the full system. The above sys- 
tem of equations is a linear system of equations for the vector of coefficients 
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The linear system can be written as 

Bo.™ = e (4) 

where e = > • • • ; dEt dt )- This system of equations can provide us 

with the time evolution of the vector 

The determination of coefficients for the reduced model through the sys- 
tem is a time-dependent version of the method of moments. We specify 
the coefficients of the reduced model so that the reduced model reproduces 
the rates of change of a finite number of moments of the solution of the orig- 
inal system. This ensures that each term in the model is properly weighted 
so that the resulting reduced model reproduces, at the scales accessible to 
the reduced model, the dynamics (see ([3])) of the original system. 

By construction, the entry By, i,j = 1, . . . , m of the matrix B measures 
the contribution of the j-th term of the reduced model to the rate of change 
of Ei. In fact, the j-th. column of the matrix B is comprised of all the 
contributions of the j-th term in the reduced model to the rates of change 
of the different Ei. While the reduced system has no need to transfer activity 
from the resolved to the unresolved scales, the columns of B corresponding 
to the activity-transferring terms will be zero (to the numerical precision 
used). Thus, the matrix B will be singular. This can be monitored by 
estimating the rank of the matrix through the Singular Value Decomposition 
(SVD) [16]. When the smallest singular value becomes non-zero for the 
numerical precision used the reduced system starts transferring activity to 
the unresolved scales. After that instant we can use the system Q to 
estimate the coefficient vector (see |14| for more details). 

The accurate determination of the coefficients guarantees that the rate 
at which the reduced model transfers activity to the unresolved scales is in 
agreement with the rate of transfer dictated by the full system. As we will 
see in Section [21 the accurate determination of the coefficients is crucial and 
moreover it can be hard. This is because the matrix B can turn out to be 
extremely ill-conditioned. Such ill-conditioning is to be expected because 
the contributions of the different terms in the reduced model to the rate 
of change of the quantities Ei can vary dramatically from term to term. 
In fact, it is this significant variation of the contributions of the various 
terms which gauges the importance of the different terms in the reduced 
model. Eventually, the difference of the importance of the different terms 
is reflected on the values of the coefficients of the terms. This will become 
apparent when we present the numerical results in Section [2j 
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1.3 The Mori-Zwanzig formalism 

We have presented in the previous section an abstract way of writing the 
reduced system which does not make any reference to a specific method for 
obtaining the functions R k (t, u'(t)) k € F appearing on the RHS of ([2]). In 

order to proceed we need to specify the functions R^\t,u'(t)). We will do 
that through the Mori-Zwanzig formalism 0[6]. 
Suppose we are given the full system 

^ =*(*,«(*)), (5) 

where u = {{u k }), k £ F U G with initial condition u(0) = no- The system 
of ordinary differential equations we are asked to solve can be transformed 
into a system of linear partial differential equations 

Lcj) k , 4> k (u ,0) = u fc, k G F U G (6) 



dt 



where L = X^eFuG ^i( u o)'^~- ^ ne solution of ([6]) is given by u k (uo,t) = 
(j>k(uo,t). Using semigroup notation we can rewrite ([6]) as 

d 

T-e tL u k = Le tL u ok 
dt 

Suppose that the vector of initial conditions can be divided as no = (uo, uq), 
where uq is the vector of the resolved variables and uq is the vector of the 
unresolved variables. Let P be an orthogonal projection on the space of 
functions of uq and Q = I — P. We delay the specification of the projection 
operator until Section [2] where we present applications to specific PDEs. 
Equation ([6]) can be rewritten as 

8 /"* 

—e tL u ok = e tL PLu ok + e tQL QLu ok + / e {t ~ s)L PLe sQL QLu ok ds, k € F, 
dt ' ' ./o 

(7) 

where we have used Dyson's formula 



e tL = e tQL + f e {ts)Lp Le sQL ds 

Jo 



(8) 



Equation ([7]) is the Mori-Zwanzig identity. Note that this relation is exact 
and is an alternative way of writing the original PDE. It is the starting 
point of our approximations. Of course, we have one such equation for each 
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of the resolved variables u k ,k € F. The first term in ([7]) is usually called 
Markovian since it depends only on the values of the variables at the current 
instant, the second is called "noise" and the third "memory". 
If we write 

e tQL QLu ok = w k , 
Wk(uo,t) satisfies the equation 

§- t w k {uQ,t) = QLw k (u , t) ^ 
w k (u , 0) = QLx k = R k (u ) - (PR k )(u ). 

If we project ([9]) we get 

d 

Pw;W k (uo,t) = PQLw k (u ,t) = 0, 
at 

since PQ = 0. Also for the initial condition 

Pw k (u ,0) = PQLu ok = 

by the same argument. Thus, the solution of (j9]) is at all times orthogonal 
to the range of P. We call ([9]) the orthogonal dynamics equation. Since 
the solutions of the orthogonal dynamics equation remain orthogonal to the 
range of P, we can project the Mori-Zwanzig equation (J7]) and find 

^-Pe tL u ok = Pe tL PLu ok + P I e {t ~ s)L PLe sQL QLu ok ds. (10) 
ot Jo 

In order to proceed with the computation of the reduced model we need 
to compute the Markovian term and the memory term. While the Marko- 
vian term is usually rather straightforward to compute (see Section [2|), the 
memory term computation is rather involved due to the presence of the or- 
thogonal dynamics evolution operator e t( ^ L . In fact, it is the presence of this 
operator which makes, in general, the computation of MZ reduced models 
prohibitively expensive (see [3] for a thorough discussion). 

One can start from (I10p and based on assumptions derive simplified 
reduced models that are easier to calculate El 13 [17]. However, it has 
proven difficult to justify the assumptions underlying the derivation of the 
simplified models. In addition, it is hard to work with the expression for 
the memory term given in (jlOp . We offer here an alternative presentation of 
the memory term which allows us to fit the MZ formalism in the framework 
described in Section [LT1 
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Note that the memory term f£ S ^ L PLe s( ^ L QLuQ k can be written, 
through Dyson's formula ©, as 

f e it - s)L PLe sQL QLu ok ds = e tL QLu ok - e tQL QLu ok (11) 
J o 

which can be rewritten (using the linearity of the operator e tL ) as 

e tL QLu ok - e tQL QLu Qk = e tL (QLu ok - e - tL e tQL QLu ok ) . (12) 
Since I = P + Q we find 

I e^ L PLe sQL QLu ok ds = e tL (QLu ok - e c ^QLu ok ) (13) 
J o 

where C(t,uo) = —tPL + [—tL,tQL] + . . . is given by the BCH series and 
[-tL,tQL] = -tLtQL - tQL(-tL) (see e.g US]). We have [—tL,tQL] = 
[tL,tPL] = [tQL,tPL]. In general, [tQL,tPL] ^ 0. However, depending on 
the initial conditions, [tQL, tPL] may be small and thus allow the simpli- 
fication of the memory term expression. The magnitude of [tQL, tPL] will 
be examined in detail in a future publication. 

If we assume that C(t, uo) ~ —tPL, expansion of the operator e~ tPL in 
Taylor series around t = gives 

P I e {t - s)L PLe sQL QLu ok ds^ (14) 
J o 

J2(-l) j+1 -Pe tL (PLyQLu ok . 

3=1 J ' 

One can obtain different simplified models by truncating the series in (|14p 
after different values of j. In particular, if we omit all the terms after the 
first one we obtain the t-model which has been studied thoroughly [61 E] - 
The series representation of the memory term in (I14p is based on the as- 
sumption of analyticity in time of the operator e~ tPL . This assumption may 
be true for small t but it does not have to hold for larger t. In other words, 
the Taylor expansion of the operator e~ tPL has, in general, only a finite ra- 
dius of convergence. Insisting on using the Taylor expansion of the operator 
e -iPL ag j g £ Qr i a £ er times is dangerous and can lead to the instability of the 
reduced model (see also Section \2A\i . In fact, when dealing with full systems 
coming from discretizations of singular PDEs, the breakdown of the Taylor 
expansion of the operator e~ tPL is related to the onset of underresolution 
on the part of the full system. 
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To proceed we need to put the MZ model given by (jlOp and ()14[) in the 
framework of Section 11.11 To do that we set 

= Pe tL PLu ok , (15) 
= (-^rfrTyP^iPLy^QLuok, j = 2,.... (16) 

With this identification we have, in essence, embedded the reduced models 
derived through the MZ formalism in a larger class of reduced models which 
share the same functional form with the MZ models but who are allowed to 
have different coefficients. In the notation of Section ll.l| the original MZ 
models correspond to the coefficient vector a^- 1 ' = (1, 1, 1, ...). 

While the original MZ models may suffer from instabilities, the new 
models can be made stable by assigning to each term in the reduced model 
the appropriate coefficient. The magnitude of the coefficient of a term re- 
flects the importance of the term in the reduced model. The values of the 
coefficients can now be determined by solving the linear algebraic system 
dH). This ensures that the coefficient of each term in the model is properly 
redefined (renormalized) so that the resulting reduced model reproduces, at 
the scales accessible to the reduced model, the dynamics (see Q) of the 
original system. 

1.4 The renormalized Mori-Zwanzig (rMZ) algorithm 

We are now in a position to state the renormalized Mori-Zwanzig algorithm 
which constructs a reduced model with the necessary coefficients computed 
on the fly. 

Renormalized Mori-Zwanzig Algorithm 

1. Choose a number of terms, say m, to keep at the Taylor expansion of 
the memory term. 

2. Evolve the full system and compute, at every step, using the SVD, the 
rank of the (m + 1) x (m + 1) matrix B. 

3. When the smallest singular value cr m+ i = TOL, (we assume that the 
singular values are indexed from largest to smallest) where TOL a 
prescribed tolerance, solve the system @ for the coefficients. 

4. For the remaining simulation time, switch from the full system to the 
reduced model with the estimated values of the coefficients. 



12 



In order to apply the algorithm above, we need to specify the quantities 
Ei, i = 1, . . . , m. Also, we need to compute the expression for the Markovian 
term, as well as the expressions for the terms in the Taylor expansion of the 
memory term. This can be a formidable task. However, as we will show in 
the next section, for the case of Burgers and Euler, these computations can 
be automated and thus we can compute arbitrarily high order terms in the 
memory term expansion with minimal effort and storage. 



2 Application of rMZ to ID Burgers and 3D Euler 
equations 

In this section we present results of the rMZ algorithm for the ID Burgers 
and the 3D Euler equations. 

2.1 ID Burgers equation 

2.1.1 Setup of the reduced model 

We use the ID inviscid Burgers equation as an instructive example for the 
constructions presented in this section. The equation is given by 

u t + uu x = 0. (17) 

Equation (|17p should be supplemented with an initial condition u(x, 0) = 
uq(x) and boundary conditions. We solve (|17p in the interval [0, 2tt] with 
periodic boundary conditions. This allows us to expand the solution in 
Fourier series 

u M (x,t) = u k (t)e ikx , 
keFuG 

where F U G = [— 4^-, 4^- — 1]. We have written the set of Fourier modes 
as the union of two sets in anticipation of the construction of the reduced 
model comprising only of the modes in F = [— y, ^ — 1], where N < M. 
The equation of motion for the Fourier mode Uk becomes 

-#■ = -2 ^ <w ( 18 ) 

p+q=k 
p,q€FUG 

To conform with the Mori-Zwanzig formalism we set 

Rk(u) = -y ^ u p u i 

p+q=k 
p,q£FUG 
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and we have 

d -^ = Rk(u) (19) 



for k G F U G. The system (|19p is supplemented by the initial condition 
uo = (uo,uq) = (uo,0). We focus on initial conditions where the unresolved 
Fourier modes are set to zero. We also define L by 

keFUG W 

Note that Lu ok = R k (uo). 

We also need to define a projection operator P. For a function h(uo) of 
all the variables, the projection operator we will use is defined by P(h(u)) = 
P(h(uQ,uo)) = /i(no,0), i.e. it replaces the value of the unresolved variables 
uq in any function h(uo) by zero. Note that this choice of projection is 
consistent with the initial conditions we have chosen. Also, we define the 
Markovian term 

R[ 1} k(u ) = PLu ok = PR k (u ) = -y ^2 uopUo q . 

p+q=k 
p,q£F 

The Markovian term has the same functional form as the RHS of the full 
system but is restricted to a sum over only the resolved modes in F. The 
full system conserves the energy | SfceFuG \ u k\ 2 contained in all the modes. 
Similarly, the Markovian term of the reduced model does not alter the energy 
content of the resolved modes. The necessary energy transfer out of the 
resolved modes rests on the memory terms. Based on our choice of projection 
operator and the scaling symmetries of the Burgers equation we set iV = 4^ . 
With the definition of P given above, we find for QLuQ k 

QLu 0k = -y U pU q - — U 0p U q - y } U pU 0q - 

p+q=k p+q=k p+q=k 

p£F peG peG 

qeG q£F qeG 

The expression for QLuQ k contains three terms which involve at least one 
wavenumber in the unresolved range G. The terms in the Taylor expansion 
of the memory term are given by 



R fk = (-lYjf^yPe'HPLy-'QLuok, j = 2, 
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For the j-th term we have 



(PL) j 1 QLu ok = (PL) j E n o P «Og-y E ^Op^Og-y E u 0p u 0q j . 

(20) 



p+q=k p+q=k p+q=k 

peF peG peG 

qeG qeF qeG 



2.1.2 Recursive computation of the memory terms 

The expression for (PLy~ 1 QLuok m (|20p for the jth term {j = 2, . . .) can be 
computed recursively using a simple construction based on a Pascal triangle. 
Note that for our choice of projection operator P, we have 

(^) i_1 (-y E ^ P u 0q )=o. 

^ p+q=k ' 
peG 
qeG 

We begin with the (first-order) term for j = 2 i.e., PLQLu^. We find 

PLQLu ok = -2— PuopPLuoq. (21) 

p+q=k 

peF 
qeG 

The first order term can be computed by convolving the resolved part of 
Puq p with the unresolved part of PLuo q . In practice, all the convolutions 
sums can be computed using Fast Fourier Transforms [19]. Note that the 
expression Puq p is linear in the Fourier modes while PLuo q is quadratic. 
Thus, the convolution sum in PLQLuok (including the factor — y)can be 
denoted as (lr * 2u), where * stands for convolution while r and u stand for 
the resolved and unresolved parts. This notation facilitates the recognition 
of the pattern for the higher order terms. With this notation, the first-order 
term can be written as 

PLQLuo k = 2 x l(lr * 2u) (22) 

where we have used bold face to denote the coefficient. We continue with 
the second order term PLPLQLu^. We find 

PLPLQLu Qk = 2(- % ^ Pu 0p PLPLu Oq -j ^ PLPu 0p PLu 0q 

p+q=k p+q=k 
peF peF 

qeG qeG 

(23) 
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The convolution sums in this term can be denoted as (lr*3u) and (2r*2u). 
The second order term can be written as 

PLPLQLuQ k = 2 x (l(lr * 3u) + l(2r * 2u)\ (24) 

To see the pattern more clearly we need one more term, PLPLPLQLuok- 
We find 

PLPLPLQLuok = (25) 
— PuopPLPLPLuoq 

p+q=k 
peF 

qeG 

-2y Yl PLPuopPLPLuoq 

p+q=k 
peF 
qeG 

~ ^ PLPLPuopPLuoq 

p+q=k 
peF 

qeG 

The terms in the parenthesis can be denoted as (lr * Au), (2r * 3u) and 
(3r * 2u). The third order term can be written as 

PLPLPLQLuok = 2 x (l(lr * 4u) + 2(2r * 3u) + l(3r * 2u) J (26) 



By examining the expressions in (I22j) - (l26p we see that the memory terms 
can be computed as weighted sums of convolution sums where the weights 
are given by appropriate Pascal triangle coefficients (the bold face numbers). 
This was to be expected since we started with a convolution sum (involving 
products) of two functions and each new term in the Taylor series involves 
a differentiation . Moreover, the number of convolution sums that need to 
be added is equal to the order of the memory term in the Taylor expansion. 
Also, for each term, the convolution sums involve expressions whose degree 
(in Fourier modes) follows an easily discernible pattern. For the Ith order 
term in the Taylor series we need the convolution sums (lr * (I + l)u), (2r * 
lu), ... ,(lr * 2u). 

Finally, the expressions entering the convolution sums can also be com- 
puted by a Pascal triangle construction. For example, in order to calcu- 
late the third order term, as can be seen from (|25p . one needs to com- 
pute and store only the quantities Puo p , PLuo p , PLPLuo p , PLPLPLuo p for 
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Figure 1: ID Burgers equation. Comparison of the evolution of the energy 
in the resolved modes computed by the random choice method, the t-model 
and various rMZ models (see text for details). 

p € FUG. Note that the quantities Puq p , PLPuq p and PLPLPuq p which are 
also needed are the same as Puo p , PLuq p , PLPLuq p for the resolved modes 
and zero for the unresolved modes. So, they do not need to be stored. They 
can be quickly constructed when needed. We see that the storage require- 
ments for the calculation of the memory terms grows only linearly in the 
order of the Taylor expansion. Also, the ability to calculate the needed ex- 
pressions through FFTs speeds up significantly the calculation of the various 
memory terms. 

The recursive estimation of the memory terms allows us to calculate 
memory terms of very high order efficiently, without having to write down 
explicitly the analytical expressions which become very complicated after 
the first few orders in the expansion. 

2.1.3 Results using rMZ 

Figure [T] shows the evolution of the energy ^ X^fcgF \ u k^ °f the resolved 
modes computed by reduced models of different orders and the random 
choice method [20]. The initial condition is Uq(x) = s'mx. This leads to 



17 



the formation of a standing shock at T = 1. All the reduced models use 
N = 16 Fourier modes while the full system has M = 32 modes. The 
results of the reduced models are compared to a converged solution of the 
random choice method with N = 4096 points. The energy of the random 
choice method solution was computed using only N = 16 modes. However, 
note that practically all the energy of the random choice method solution 
is concentrated in the first few Fourier modes, so even if we had computed 
the energy for all N = 4096 Fourier methods the results would not have 
changed. This is to be expected, since for the initial condition we are using, 
a standing shock forms at time T = 1 and, thus, by time T = 100 the only 
Fourier modes having some energy left in them are the first few. 

The quantities E{ used to set up the linear algebraic system needed to 
compute the coefficients of the reduced system are l p norms of the solution. 
In particular, for the first order model we use Ei = YlkeF l M fc| 2 *) i = 1, 2. For 
the second order model Ei = YLk&F l M fc| 2 \ * = 1, 2, 3 and for the third order 
model Ei = X^fceF l n fc| 2 *> * = 1,2,3,4. In general, for the reduced model 
of order A we need A + 1 quantities because we also have to compute the 
coefficient of the Markovian term. All the calculations are done in double 
precision. The tolerance TOL which is used to decide when it is time to 
switch to the reduced model is set to TOL = 10~ 12 . The systems of ordinary 
differential equations for the different reduced models were solved using the 
Runge-Kutta-Fehlberg method with the stepsize control tolerance set to 
10" 10 p]. 

The numerical problem of solving the linear system for the coefficients is 
hard because the resulting system has very large condition number and very 
small determinant. This happens for three reasons. First, the dominant 
contribution to the linear system matrix comes from the Markovian term 
(except for the contribution to the rate of change of the E\ = YlkeF \ u k\ 2 
which is zero). This means that the coefficient of the Markovian term is 
practically 1. Second, the contributions of each memory term to the rates of 
change of the different Ei vary dramatically. Third, the contributions to each 
Ei by the different memory terms also varies substantially. Of course, this 
situation is exacerbated if we use more terms we use in the expansion. For 
the case when we retain up to the third order term in the memory expansion, 
we have to deal with condition numbers of the order 10 11 and determinant 
values of order 10~ 20 . Inevitably, even the use of double precision cannot 
provide us with an accurate estimate of the coefficients. Since the linear 
system matrix is practically singular (for the numerical precision used) we 
have chosen to solve the linear system using the SVD algorithm |16j . 
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A partial remedy to the problem comes from a slight modification in 
the way of estimating the coefficients. Since we know that the Markovian 
term coefficient is practically 1, we can set it to 1, and subtract the column 
of contributions of the Markovian term from the RHS of the linear system. 
This allows us to reduce the dimensionality of the linear system to be solved 
from (A + 1) x (A + 1) to A x A. This practice of subtracting almost equal 
numbers is not advisable in general because it leads to loss of significant 
digits |16j . However, in our case it helps to improve the results by lowering 
the condition number of the matrix from about 10 11 to about 10 5 . In Figure 
1, the estimation of the coefficients for the third order model using the 
reduced dimension matrix is denoted as "rMZ-reduced 3rd". 

As shown in Figure Q3 the rMZ models of first and second order give 
practically the same results as the t-model. The third order model gives a 
slight improvement. However, when the reduced dimension matrix is used, 
the energy evolution predicted by the third order model is practically iden- 
tical to the correct energy evolution of the resolved modes predicted by the 
random choice method. If we increase the resolution of the reduced model, 
the numerical problems for the calculation of higher order coefficients be- 
come even more pronounced. This is to be expected, since a larger resolu- 
tion means that the renormalized coefficients of the reduced model will be 
smaller. Thus, computing them with accuracy is more difficult. 

We have to note, that the computation of the higher order coefficients 
is more difficult for our choice of initial condition since it only involves one 
active Fourier mode. Initial conditions with more active Fourier modes will 
transfer activity to the unresolved scales at a higher rate and thus the cor- 
responding renormalized coefficients will be larger. A detailed study of the 
behavior of the coefficients for different initial conditions will be presented 
elsewhere. For the first order model, we have already presented in p3] a 
detailed study about the change of the value of the renormalized coefficient 
with resolution up to the order of 10 5 Fourier modes. In that work, the 
renormalized coefficient calculation was used to determine, in a fixed point 
analysis, the blow-up exponent. 

2.2 3D incompressible Euler equations 

Consider the 3D incompressible Euler equations with periodic boundary 
conditions in the cube [0,27r] 3 : 

u t + u-Vu = -Vp, V • u = 0, (27) 
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where u(x,t) = (ux(xi, x 2 , x 3 , t), u 2 (xi, x 2 , x 3 , t), u 3 (x\, x 2 , x 3 , t)) is the ve- 
locity, p is the pressure and V = (gfj, gfj, gf^)- The system in (|2.2p is sup- 
plemented with the initial condition u(x,0) = uq(x) which is also periodic 
and incompressible and x = (x\,x 2 ,x 3 ). Since we are working with periodic 
boundary conditions, we expand the solution in Fourier series keeping N 
modes in each spatial direction, 



u M (x,t) = u k (t)e 



ikx 

fceFuG 



where FUG = [-f ,f -l]x[-f ,f -l]x[-f ,f -1]. Also A; = (k u k 2 ,k 3 ) 
and u k (t) = (ul(t),ul(t),u 3 k (t)). 

The equation of motion for the Fourier mode becomes 

= -i ^ k ' U P A kU q , (28) 
p+q=k 
p,q£FUG 

where Aj, = I — jj^j is the incompressibility projection matrix and I is the 
3x3 identity matrix. The symbol • denotes inner product in M 3 . The system 



([28]) is supplemented by the initial condition uo = {«&(())} = {uq^}, k € 
FUG, where UQk are the Fourier coefficients of the initial condition uq{x). 
To conform with the MZ formalism we set 



p+q=k 
p,qtEFUG 



= (29) 



and we have 



for k € F U G. The system (J29J) is supplemented by the initial condition 
uo = (uq,uq) = (no,0). Note that we focus on initial conditions where the 
unresolved Fourier modes are set to zero. We also define L° by 



d 

L = ^2 R k(uo) 



dunk 

keFUG u/c 

Note that Luofc = -Rfc(^o)- Consider the subset F = [— 1] x [— y, y — 

1] x [— y, y — 1] for N < M. We will construct the reduced models for the 
Fourier modes u k with k S F. 
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We need to define a projection operator P. For a function h(uo) of all 
the variables, the projection operator we will use is defined by P(h(u)) = 
P(h(uQ,uo)) = h(uo,0), i.e. it replaces the value of the unresolved variables 
no in any function h(uo) by zero. Note, that this choice of projection is 
consistent with the initial conditions we have chosen. Based on our choice 
of projection operator and the scaling symmetries of the Euler equations we 
set N = f. 

Define 

Rk(uo) = PRk(uo) = ~i ^ k ' UOpAkUOq- 

p+q=k 
p,q€F 

The Markovian term has the same functional form as the RHS of the full 
system but is restricted to a sum over only the resolved modes in F. The 
full system conserves the energy \ J2keFuG \ u k\ 2 contained in all the modes. 
Similarly, the Markovian term of the reduced model does not alter the energy 
content of the resolved modes. The necessary energy transfer out of the 
resolved modes rests on the memory terms. 

With the definition of P given above, we find for QLuo k 

QLu ok = -i ^2 k- u pA k u q - ^ A; ■ u pA k u 0q - i ^ k ■ uo p A k u 0q . 

p+q=k p+q=k p+q=k 

p&F peG peG 

q£G q£F qt=G 

The expression for QLuo k contains three terms which involve at least one 
wavenumber in the unresolved range G. The terms in the Taylor expansion 
of the memory term are given by 

R fk = (-l) 3 ^^yPe tL {PLy- l QLu ok , j = 2,... 

For the j-th term we have 

(PL) j ~ 1 QLu ok = (PL) J_1 f —i y~] k ■ u 0p A k u 0q - i ^ k ■ u 0p A k u 0q 

p+q=k p+q=k 
peF p€G 
qeG qGF 

(30) 



^ k ■ uo p A k u 0q j . (31) 



p+q=k 

P eG 
q eG 
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Figure 2: 3D Euler equation. Comparison of the evolution of the energy in 
the resolved modes computed by the t-model and various rMZ models (see 
text for details). 

The different terms in the memory expansion can be computed recursively 
as in the case of Burgers. However, there is a slight complication because 
the presence of the incompressibility operator and of the inner product on 
the RHS destroys the commutativity which allowed us in Burgers to group 
terms (the factor 2 which appears outside every parenthesis there). This 
problem can be addressed by a construction which uses 2 Pascal triangles 
instead of 1 used in the case of Burgers. For each order, one adds up the 
corresponding terms from the 2 Pascal triangles and obtains the desired 
memory term. Other than that, the recursive algorithm remains the same 
and we omit the details. Also, note that all the higher order terms are 
divergence-free by construction. 

We have used the same quantities Ei as in the case of Burgers, with the 
obvious generalizations, since for 3D Euler we have a 3-dimensional velocity 
vector instead of the scalar velocity in Burgers. Also, even though for 3D 
Euler we have a 3-dimensional vector, we have assumed that the reduced 
model renormalized coefficients are the same for all 3 velocities. This is 
a simplifying assumption. Of course, one can use different renormalized 
coefficients for the different velocities at the expense of having to solve a 
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larger linear system for the renormalized coefficients. A detailed study of 
that case will be presented in future work. 

We have used the Taylor-Green initial condition (see e.g. [22] and refer- 
ences therein) which is given by 

ui(x,0) = sin(xi) cos(x2) cos(x3), 
U2(x,0) = — cos(xi) sin(x2) cos(x3), 
u 3 (x,0) = 0. 

Figure [2] shows the evolution of the energy | J2keF \ u k\ 2 of the resolved 
modes for different resolutions computed by rMZ reduced models of different 
orders and the t-model. We have presented results for the rMZ reduced 
models using the reduced linear system matrix approach discussed above to 
tame the condition number of the matrix. Based on these results, we make 
two observations. 

First, for 8 3 resolved modes, the rMZ third order model dissipates energy 
at a slower rate than the i-model with 8 3 resolved modes. This is true not 
only for the third order model but also for the first and second order models 
(we have omitted those results to avoid cluttering the figure). This slower 
rate of energy dissipation compared to the i-model holds also for the case 
of 16 3 resolved modes. 

The second observation is that the rate of energy dissipation of the rMZ 
models is consistent with the rate predicted by the t-model with higher 
resolution. This is to be expected since a higher order model should result 
in a more accurate prediction of the energy dissipation rate. 

The reader may be concerned about the small resolutions used in the 
numerical experiments. There are two reasons for that. First, if one keeps 
several terms in the memory expansion, then, for a very smooth initial con- 
dition like the one we use, the matrix B becomes even more ill-conditioned 
for large resolutions. However, this is not a severe problem. On the contrary, 
it signifies that most of the higher order terms should have small coefficients 
and thus can be safely removed from the model. 

The second reason we have used small resolutions both for Burgers and 
Euler is because an accurate reduced model should be able to reproduce the 
correct energy content for its resolved scales no matter how small the reso- 
lution. For example, for Burgers, where we know what the energy content 
should be after the singularity, we see that the rMZ model with a small reso- 
lution (16 3 ) indeed reproduces the correct energy content for this resolution. 
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Ratio of smallest active scale in initial condition to resolution of reduced model 



Figure 3: Comparison, for ID Burgers and 3D Euler, of the value of the 
renormalized coefficient for the first order rMZ model for different resolu- 
tions. The initial conditions for Burgers and Euler have only 1 active Fourier 
mode in each direction (see text for details). 

2.3 Universality of the renormalized coefficients 

We have hinted in the introduction at the possibility that the renormalized 
coefficients may be determined by two factors: i) the ratio of the smallest 
active scale in the initial condition to the smallest resolvable scale and ii) 
the scaling symmetries of the equation under investigation. Even though 
the numerical difficulties with the ill-conditioned linear system matrix do 
not allow us at present to study accurately the higher order renormalized 
coefficients, we have enough accuracy to study the first order renormalized 
coefficient both for ID Burgers and 3D Euler. Note that the two equations 
share the same scaling symmetries. Also, we have chosen for Burgers the 
initial condition uo(x) = sinx which has only one active Fourier mode (for 
k = ±1) and for 3D Euler, the Taylor-Green initial condition which also has 
only active Fourier modes for fej = ±1, i= 1,2, 3. 

Figure [3] shows the comparison of the value of the renormalized first-order 
coefficient for Burgers and 3D Euler as a function of the ratio of the smallest 
scale active in the the initial condition to the smallest scale of the reduced 
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model. We make two observations. First, the values of the renormalized 
coefficient for the two equations are in remarkable agreement. Second, from 
the slope of the linear fit, we see that the value of the coefficient is practically 
equal to the ratio of the smallest scale active in the the initial condition to 
the smallest scale of the reduced model. 

Needless to say that one example is not enough to infer the generality of 
the result for arbitrary initial conditions. A theoretical explanation of this 
result is lacking at the moment. Note that due to the way we have defined 
the terms in the expansion of the memory, all the terms have the same di- 
mensions as the Markovian term and the left-hand side of the equation for 
each Fourier mode. So, the corresponding coefficients have to be dimen- 
sionless. Thus, we expect the coefficients to depend on ratios of quantities 
with the same dimensions. Here we have investigated the possibility that 
this ratio is that of the smallest active scale in the initial condition to the 
smallest active scale of the reduced model. 

We should comment here on the behavior of the rMZ algorithm for the 
2D Euler equations for which the 2D version of the Taylor-Green initial 
condition is an exact solution, i.e. a steady state. Exactly because it is a 
steady state there is no need for a reduced model. Application of the rMZ 
algorithm agrees with this. There is never any need to transfer energy to 
the unresolved scales and thus, no need to switch to a reduced model. The 
contributions of the different memory terms to the matrix B are all well 
below the double precision threshold. This allows the freedom to assign to 
the renormalized coefficient the values shown in Figure [3] without incurring 
any trouble. In other words, the behavior of the solution of the 2D Euler for 
the Taylor-Green initial condition does not contradict the agreement for the 
renormalized coefficient of the ID Burgers and 3D Euler equations shown in 
Figure [3l 

2.4 rMZ vs MZ 

In this section we present results which show that the renormalized version 
of the MZ formalism is advantageous with respect to the original MZ formal- 
ism. In particular, we show that for the same order in the Taylor expansion 
of the memory term, the renormalized algorithm leads to the stabilization 
of the reduced model. Figure H] compares, for the ID Burgers equation, 
the energy i X^fceF l^fcl 2 f° r 16 resolved modes for the renormalized and 
unrenormalized third order models. Figure [5] compares, for the 3D Euler 
equations, the energy ^ SfceF \ u k\ 2 for 8 3 resolved modes for the renormal- 
ized and unrenormalized third order models. The unrenormalized models 
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Figure 4: Comparison, for ID Burgers, of the value of the energy content 
of the resolved modes for the third order renormalized and unrenormalized 
MZ models. 

quickly become unstable and lose all predictive ability. 

The unrenormalized expansion leads to divergence of the predicted en- 
ergy content of the resolved modes. This is analogous to the divergences that 
plagued perturbative calculations in quantum field theory (QFT) before the 
advent of renormalization [23] . In QFT, the reason for the divergences was 
that the perturbation expansion was performed in a quantity (the bare mass, 
charge etc.) which turns out to be ill-defined. The process of renormaliza- 
tion replaces the perturbation expansion in powers of the ill-defined quantity 
with a perturbation expansion in powers of the experimentally determined 
values of this quantity. This allows the subtraction of the terms that cause 
divergences and leads to finite results. 

In our case, the expansion of the memory of the MZ formalism is ill- 
defined because the Taylor expansion at t = breaks down after some time. 
On the other hand, the renormalized MZ formalism takes into account dy- 
namic information from the evolution of the full system (while this system 
is still valid) and prescribes to each term in the memory expansion an ap- 
propriate coefficient. The coefficient measures how important this term is. 
In this way, the divergences are averted and the results become finite. 



26 




Figure 5: Comparison, for 3D Euler, of the value of the energy content of 
the resolved modes for the third order renormalized and unrenormalized MZ 
models. 
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3 Conclusions and future work 



We have presented a new way of computing reduced models for systems 
of ordinary differential equations. The approach combines renormalization 
and effective field theory techniques with the Mori-Zwanzig formalism. The 
constructed reduced models are stable because they transfer activity out 
of the resolved scales at a rate which is dictated by the full system. The 
consistency between the rate of transfer activity of the reduced model and 
the rate of transfer activity dictated by the full system is the analog of 
the matching conditions employed in effective field theory. The matching 
conditions lead to a redefinition (renormalization) of the coefficients of a 
reduced model originally constructed through the Mori-Zwanzig formalism. 

The results we have obtained for the ID Burgers and 3D Euler equations 
are rather encouraging. The proposed approach of calculating the renormal- 
ized coefficients leads to a linear algebraic system for the coefficients. The 
matrix of this linear system is ill-conditioned by construction. The reason 
is that the reduced model comprises of two parts. The Markovian part 
which is of the same functional form as the full system but defined only 
on the resolved scales, and the memory part which has no analog in the 
full system. The memory part is there to compensate for whatever activ- 
ity the Markovian term cannot reproduce. Thus, the determination of the 
coefficients appearing in the memory part are, in essence, computed by sub- 
tracting the contribution of the Markovian term from the full system. This 
difference may be small and also the different terms in the memory part 
can have varying contributions to the difference. These factors result in the 
ill-conditioning of the linear system that needs to be solved. In the future, 
we plan to address the problem through various techniques designed to deal 
with ill-conditioned matrices. 

The results for ID Burgers and 3D Euler equations also point to the pos- 
sibility that the renormalized coefficients may be completely determined by 
the structure of the initial condition and the scaling symmetries of the equa- 
tion. This needs to be investigated further and a detailed study involving 
random initial conditions with various spectra will be presented elsewhere. 

The proposed approach can also be applied to the Navier-Stokes equa- 
tions [JJ. It is straightforward to incorporate viscosity in the new expansion 
of the memory term proposed in the current work. We note that the vis- 
cosity starts contributing from the second order memory term. Also, the 
inclusion of viscosity does not complicate considerably the recursive algo- 
rithm for the calculation of the higher order terms. The expressions needed 
to compute the viscosity contributions can be estimated through terms al- 
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ready computed in the Pascal triangle construction for the inviscid terms. 

The approach presented in the current work opens new possibilities for 
the construction of accurate and stable reduced models for (large) systems 
of ordinary differential equations. It also highlights the affinity between 
problems of model reduction in scientific computing and the construction of 
effective field theories in high energy physics. We hope that this connection 
will benefit the problem of constructing reduced models and will be of use in 
tackling real world problems which are impossible to address through brute 
force calculations. 

In conclusion, as Steven Weinberg once put it [21], renormalization is 
indeed a good thing. 
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